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ABSTRACT 


Various  methods  for  numerically  generating 
orthogonal  boundary-fitted  coordinate  systems  are  dis¬ 
cussed.  Simple  generating  equations  are  used  with 
coordinate  distributions  specified  on  all  of  the  boun¬ 
daries.  Iterative  determination  of  orthogonalizing 
source  terms  for  Poisson  generating  equations  con¬ 
verges  very  slowly.  An  iterative  solution  of  a  new 
generating  system  yields  better  results  and  is  applied 
successfully  to  several  geometries. 


ADMINISTRATIVE  INFORMATION 

This  work  was  performed  under  NAVSEA  Mathematical  Sciences  Program,  "Numerical 
Methods  for  Naval  Vehicles,"  Program  Element  61153N,  Task  Area  SR0140301,  Task 
15321,  DTNSRDC  Work  Unit  1808-010. 


INTRODUCTION 

Conformal  mappings  which  are  analytic  functions  of  a  complex  variable  have 
found  extensive  application  to  problems  in  physics.  However,  these  mappings  are 
restricted  to  two  dimensions  and  have  other  limitations  which  sometimes  seriously 
diminish  their  usefulness. 

In  recent  years  more  general  transformations  have  come  into  use  under  the 

1  2* 

names  of  "boundary-fitted"  or  "surface-oriented"  coordinate  systems.  ’  These 
coordinate  systems  are  computed  as  approximate  numerical  solutions  to  elliptic  gen¬ 
erating  equations.  Usually  the  coordinate  distribution  on  the  boundaries  of  the 
physical  regions  can  be  specified  arbitrarily.  The  standard  conformal  system  is 
obtained  only  when  the  Laplace  equation  is  used  with  a  particular  boundary  coordi¬ 
nate  distribution.  Thus  boundary-fitted  coordinates  encompass  a  much  larger  set  of 
systems  (including  three-dimensional  ones)  than  the  subset  of  conformal  mappings. 
However,  these  boundary-fitted  systems  are  not,  in  general,  orthogonal.  This 

*A  complete  listing  of  references  is  given  on  page  29. 
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nonorthogonality  can  have  a  deleterious  effect  on  accuracy,  stability,  and  computa¬ 


tional  complexity. 

Some  effort  has  been  devoted  to  finding  an  intermediate  set  of  coordinate  sys¬ 
tems  which  retain  most  of  the  flexibility  of  the  general  boundary-fitted  systems 

3  4 

but  yet  are  orthogonal.  Potter  and  Tuttle  and  Ghia,  Hodge,  and  Hankey 
presented  a  method  for  orthogonalizing  a  nonorthogonal  boundary-fitted  system. 

Other  methods  for  generating  orthogonal  systems  were  presented  in  a  note  by  Mobley 

and  Stewart^  and  in  work  reviewed  by  Eisemann.^  However,  in  all  cases  the  specifi¬ 
cation  of  coordinate  distribution  was  somehow  restricted  on  at  least  one  boundary. 
The  present  report  discusses  various  methods  for  creating  orthogonal  coordinate 
systems  using  simple  generating  equations.  These  methods  allow  coordinate  distri¬ 
butions  to  be  specified  arbitrarily  on  all  of  the  boundaries.  Iterative  determina¬ 
tion  of  orthogonalizing  source  terms  for  Poisson  generating  equations  is  attempted 
with  very  limited  success.  The  iterative  solution  of  a  new  generating  system 
yields  better  results  and  successful  application  to  several  geometries  is 
presented . 


METHODS  FOR  CONSTRUCTING  ORTHOGONAL  SYSTEMS 
Many  investigators  have  used  the  system  of  Poisson  equations 

£  +  £  =  P 

xx  yy 

n  +  n  =  Q 

xx  yy 


(1) 


to  generate  numerical  transformations  from  Cartesian  coordinates  (x,y)  to  the 
transformed  coordinates  (£,n)  for  a  variety  of  geometries.  The  forcing  functions  P 

and  Q  in  Equation  (1)  provide  a  means  for  influencing  or  controlling  the  coordinate 

„  , _  ,  ,  1,2,4, 7-12 

system  obtained. 

For  ease  in  the  numerical  solution  of  Equation  (1)  all  the  calculations  are 


! 


I 
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done  in  the  transformed  plane  on  a  uniform  square  mesh.  Interchanging  the  depen¬ 
dent  and  independent  variables  gives 

axrr  -  28x.  +  yx  +  J' (Px^  +  Qx  )  =0 

££  Cn  nn  £  n 

(2 

ayCS  '  26yCn  +  Yynn  +  j2(PyC  +  QV  =  0 


where 


a  .  x2  +  y2  s  -  *txn  +  ytyn 

T  *  *E  +  A  3  "  Vn  “  Ve 


Replacing  Equation  (2)  with  central  difference  formulae  yields 


i.  . 

/  3  =  [(a.  .  +  J?  .  P 
'n  1  J  i  J  i  J 


i+l,j 


+  (a.  .  -  J?  .  P.  ./2) 

i  .1  ill] 


i_1,j  +  (y .  .  +  J?  .  Q.  . /2)  i,j+1 

'i-l.j  13  1  1  1  J 


+  (y,  .  -  J?  .  Q.  . / 2) 
i  J  i  J  i  J 


i.l-l 


-  <6i  j/2> 


Xi+l,j+l  +  Xi-1, j-1 
yi+l,j+l  +  yi-l,j-l 


i-l’J+l  t+1^~1  ]/2(«.  .  +  y.  .) 

yi-l,j+l  "  yi+l,j-l  1  J  13 


where  a  ,  6^  ,  and  are  central  difference  approximations  to  the 

corresponding  variables  of  Equation  (3). 

Extensive  use  has  been  made  of  the  Poisson  and  Laplace  (P=Q=0)  generating  sys- 

12  4  7 

terns  of  Equation  (1)  in  the  solution  of  fluid  flow  problems.  ’  ’  ’  Several 
8—10 

investigators  used  simple  exponential  forcing  functions  to  attract  lines  into 
or  repel  lines  from  certain  regions  of  the  physical  plane.  Thompson^  and 
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Plant  "  have  suggested  a  method  for  an  a  priori  determination  of  P  and  Q  which 
will  retain  boundary  spacing  throughout  the  region. 

In  this  report  we  discuss  automatic  methods  for  finding  functions  P  and  Q 
which  assure  an  orthogonal  mesh  for  a  given  geometry.  The  condition  for  ortho¬ 
gonality,  "  =  constant  lines  perpendicular  to  n  =  constant  lines,  is  $  =  0,  since 

6  =  0  =>  xr/y_  =  -v  /x 

£  £  n  n 

which  is  equivalent  to 


i/v 


1  n=constant 


-y. 


£=constant 


That  is,  the  slopes  of  the  two  sets  of  coordinate  lines  are  negative  reciprocals  of 
each  other.  An  orthogonal  coordinate  system  generated  in  this  manner  would  be  a 
solution  of  Equation  (2)  without  the  cross  derivative  terms,  i.e., 

axrr  +  yx  +  J2(Px_  +  Qx  )  =0 

££  nn  £  n  ^ 

ay^  +  Yynn  +  J2(Py.  +  Qy^)  =  0 

Two  important  questions  arise:  Under  what  conditions  do  orthogonalizing 
source  terms  P  and  Q  exist?  If  they  exist,  are  they  unique?  Since  we  do  not  now 
have  the  answers  to  these  questions,  it  is  very  important  to  have  a  test  problem 
for  which  an  orthogonal  grid  is  known  to  exist.  Failure  to  generate  an  orthogonal 
grid  for  such  a  problem  can  then  safely  be  attributed  to  the  method  used. 

A  useful  test  problem  is  illustrated  in  Figures  1  and  2.  These  figures  show  a 
rectangular  region  with  uniform  boundary  coordinate  spacing  in  the  horizontal 
direction  and  nonuniform  boundary  coordinate  spacing  in  the  vertical  direction. 

The  grid  shown  in  Figure  1  was  obtained  by  solving  Equation  (4)  with  P  =  Q  =  0  using 
successive  overrelaxation  (SOR).  Note  that  one  characteristic  of  a  Laplace- 
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generated  mesh  is  that  the  coordinate  lines  tend  toward  equal  spacing  in  the  inte¬ 
rior  for  any  boundary  spacing.  This  effect  can  be  seen  in  Figure  1  near  the  center 
of  the  region. 

An  orthogonal  grid  for  this  boundary  spacing  can  be  generated  using  Equation 
12)  with  source  terms  determined  by  the  boundary  coordinate  distribution  as  fol¬ 
lows  : 


P.  . 
13 


(-x„,/x3) .  ,  and  Q..  =  (-v  /y3),  . 

kK  t,  i,l  ij  -  nn  n  l,j 


(6) 


where  i  =  1  and  j  =  1  correspond  to  vertical  and  horizontal  boundaries.  This  form 

1? 

for  the  forcing  functions  is  derived  by  Plant  “  following  a  method  suggested  by 
11  4 

Thompson  and  Ghia.  The  orthogonal  grid  thus  generated  is  shown  after  40  SOR 
iterations  in  Figure  2.  (In  this  example  P  =  0  and  Q  =  const.  0.)  We  have 

tested  methods  for  generating  orthogonal  grids  for  their  ability  to  reproduce  the 
coordinate  system  of  Figure  2. 

One  method  of  generating  a  transformation  which  yields  an  orthogonal  coordi¬ 
nate  system  is  to  solve  the  systems  of  Equations  (2)  and  (5)  simultaneously.  These 
systems  represent  four  equations  in  the  four  unkowns ,  x,  y,  P,  and  Q.  A  solution, 
if  it  exists,  would  satisfy  =  0. 

Unfortunately,  solving  the  combined  system  of  Equations  (2)  and  (5)  is  not 
easy,  but  iterative  methods  are  applicable.  For  instance,  the  finite-difference 
versions  of  Equation  (2)  can  be  used  to  obtain  new  estimates  of  x  and  y  using 
current  values  of  P  and  Q.  Then  Equation  (5)  in  finite-difference  form  can  be  used 
to  find  new  estimates  of  P  and  Q  using  the  updated  values  of  x  and  y.  Solving 


Equation  (5)  for  P  and  Q  yields 


Pn  -  [a(x  y_.  -  v  x  3  +  y (x  y  -  y  x  )]/J3 
n  U,  '  n  fit  n  nn  qn 

<)n  =  [a(yrxrf.  -  x.yrr)  +  YCy^  -  x^)]/.]3 


A 


where  the  superscripts  indicate  iteration  level.  Solving  Equation  (2)  for  P  and 
Qn  1  yields 

Pn  1  =  [a (x  v  -  y  x  )-26(x  y  -  y  x  )+y(x  y  -  y  x  )]/J3 

n  ££  n  n Jn  n  nn  n  nn 

(8) 

Q""1  -  [«(yE*u  -  xcy5c)-2B(yE%  -  Vcn>'H<lr!;xnn  -  x^l/J5 

Comparison  of  Equations  (7)  and  (8)  shows  that 

Pn  =  p"'1  +  26(xny^  -  ynX^)/J3  (9a) 

qn  =  Qn_1  +  20(y5xCn  -  vW/j3  (9b) 

Thus,  solving  Equations  (2)  and  (5)  iteratively  is  equivalent  to  using  Equations 

(2)  and  (9)  in  an  iterative  scheme.  Note  that,  according  to  Equation  (9),  if 
n  n-1  n  n-1 

=  0,  then  P  =  P  and  Q  =  Q  and  the  iteration  process  for  P  and  Q  will  have 
converged . 

Unfortunately,  the  system  formed  by  Equations  (2)  and  (9)  cannot  handle  the 
test  problem.  This  can  be  demonstrated  by  considering  what  happens  in  an  attempt 
to  achieve  the  orthogonal  coordinate  system  of  Figure  2  by  iterating  according  to 
Equations  (2)  and  (9),  with  the  coordinate  system  of  Figure  1  as  a  starting  point. 
Consider  the  local  grid  configurations  in  Figures  3(a)  and  3(b),  both  of  which 
appear  in  Figure  1.  Equation  (9)  predicts  an  initial  Qb  (Q°  =  0)  of  different  sign 
for  the  two  configurations,  even  though  the  ultimate  Q  desired  is,  according  to 
Equation  (6),  uniform  over  the  region.  Thus  iterating  according  to  Equation  (9) 
leads  to  deterioration  of  the  system  in  part  of  the  region,  and  the  calculations 
diverge.  This  failure  is  not  of  the  basic  system  of  Equations  (2)  and  (5)  but 
rather  of  the  iterative  method  of  solving  the  system. 


Further  study  of  Figure  3  can  lead  to  an  improved  method.  Equation  (9) 
represents  an  adjustment  of  the  source  term  Q^,  depending  on  ^et  the 
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configurations  in  Figure  3  do  not,  by  themselves,  yield  any  information  about  the 
desired  movements  of  the  grid  point  (i,  j).  Rather,  Figure  3(a)  suggests  a  move¬ 
ment  of  the  point  (i-l,j)  downward  and  (i+l,j)  upward  to  improve  orthogonality. 
This  leads  to  a  modification  of  Equation  (9b)  of  the  form 


Qn+!  . 
1-1.1 

+  R  1 

IJ 

(10a) 

on+1 
•i+l,  i 

-  6  | 

IJ 

Vtn  -  Vtn)/j3|i+l,j 

(10b) 

which  will  move  points  in  a  desirable  manner  in  both  Figures  3(a)  and  3(b). 
Replacing  i  by  i+1  in  Equation  (10a)  and  by  i-1  in  Equation  (10b)  and  then  adding 
the  results  yields 


n+l  nn  +  6l+1’i  ~  gl~1»3 
Qij  '  \i  2 


<yt*Cn  '  Vtn)/j3l 


ij 


or 


(11a) 


Equation  (11a)  seems  to  provide  a  workable  scheme  along  with 


Pn.  +  R  I (x  yr  -  y  x  )/J; 
ij  V  n  Cn  n  Co 


i] 


(lib) 


which  is  arrived  at  through  similar  arguments. 

The  coordinate  system  shown  in  Figure  4,  which  is  almost  identical  to  that  of 
Figure  2,  was  generated  from  the  initial  system  of  Figure  1  in  1500  iterations 
using  Equations  (2)  and  (11).  Overrelaxation  was  used  in  the  source  term  iteration 
in  an  unsuccessful  attempt  to  accelerate  the  convergence. 

A  nonorthogonal  mesh  with  nonuniform  boundary  coordinate  spacing  in  both  vert¬ 
ical  and  horizontal  directions,  as  shown  in  Figure  5,  was  used  as  the  initial  guess 
in  a  further  test  of  the  system  of  Equations  (2)  and  (11).  In  this  case  both  P  and 
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Figure  5  -  Nonorthogonal  Coordinate  System  with  Unequal  Spacing  on  all  Boundaries 


Q  source  terms  had  to  be  determined.  The  result  after  1000  iterations  is  shown  in 


Figure  6.  The  "waviness"  of  the  mesh,  even  after  so  many  iterations,  indicates  the 
convergence  difficulties  of  this  method  which  may  preclude  it  in  many  cases  as  a 
useful  scheme  for  generating  orthogonal  grids. 

In  applying  Equations  (2)  and  (11)  we  are  attempting  to  find  a  coordinate  sys¬ 
tem  with  source  terms  P  and  Q  such  that 

6^  =  Bn  =  0  (!2) 


The  solution  0  =  constant  to  Equation  (12)  exists  only  if  it  is  consistent  with  the 
boundary  data.  Only  at  the  corners  of  the  computational  region  is  0  specified  in 
advance.  Thus  we  can  hope  for  a  solution  if  Q  is  the  same  at  all  corners  and  for 
an  orthogonal  solution  only  when  0  =  0  at  the  corners. 

A  different  approach  to  this  orthogonal  grid  problem  is  to  bypass  the  Poisson 
equations  altogether  and  to  choose  a  generating  system  based  entirely  on  0  . 

Once  the  grid  points  have  been  suitably  positioned  in  the  physical  region,  the 
forcing  functions  for  the  Poisson  system  that  would  generate  this  configuration 
could  be  found  from  Equation  (7). 

Even  though  0  =  0  is  the  condition  for  or thogonality ,  this  equation  alone  is 
not  sufficient  for  obtaining  the  required  transformation.  Two  equations  are 
needed,  since  we  must  find  both  the  x  and  y  coordinates  of  the  transformed  points. 

As  a  new  generating  system,  consider  Equation  (12).  Expanding  Equation  (12) 


gives 


xrx_  +  x_._x  +  y_y.  +  y_,.y  =  0 

£  £n  ££  n  £J£n  £f/n 


xrx  +  xr  x  +  y  v  +  y.  y  =  0 

£  nn  £n  n  £  nn  £n  n 


(13a) 

(13b) 


To  compute  the  transformation,  Equations  (13a)  and  (13b)  are  combined  as  follows: 


the  product  of  Equation  (13a)  and  x^  is  added  to  the  product  of  Equation  (lib) 
and  x yielding 


Xn  xrr  +  xr  x  +  2x.x  x  +  x  y  y__  +  x  y  y 

n  C  nn  t,  n  irn  ££  f/f/nn 

+  (x  y  +  X  y  )  y ,  =  0 

n  f,  i,  n  f,n 


(14a) 


the  product  of  Equation  (13a)  and  y  is  added  tc  the  product  of  Equation  (lib) 


and  yr,  yielding 


;n  y« 


+  v 

yS  ynn 


+  x  y  x 

n  n  f.f. 


+  x  y  x 

f,  r  nn 


(14b) 


+  (x  y  +  x  y  )x_  =  0 

n  r  f,  n  bn 

The  reason  for  replacing  Equations  (13)  with  Equations  (14)  is  to  obtain  a  non-zero 
coefficient  for  x„  and  y^  in  the  finite-difference  forms  of  Equations  (14a)  and 
(14b).  This  eliminates  the  possibility  of  dividing  by  zero  in  the  iteration  pro¬ 
cess.  Each  derivative  in  Equations  (I4a)  and  (14b)  is  replaced  by  the  appropriate 
central  difference  formula  and  the  system  is  solved  iteratively  using  successive 
overrelaxation . 

The  system  of  Equation  (12)  not  only  successfully  generates  the  orthogonal 
grid  of  Figure  2  but  also  generates  the  acceptable  coordinate  system  of  Figure  7 
from  the  initial  system  of  Figure  5  in  500  iterations,  half  the  number  of  itera¬ 
tions  needed  to  generate  the  unacceptable  system  of  Figure  6  with  Equations  (2)  and 


(11). 


FURTHER  EXAMPLES 

in  each  of  the  final  six  examples,  we  show  a  nonor thogonal  coordinate  system 
generated  by  Equation  (2)  with  P  2  Q  2  0  (a  Laplace  system)  and  a  coordinate  system 
generated  by  Equation  (12).  The  first  of  these  examples,  shown  in  Figure  8,  is  a 
simply-connected  region  with  one  convex  boundary.  Next  we  have  a  similar  region 
with  a  concave  rather  than  convex  curved  boundary  as  seen  in  Figure  9.  Note  that 
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Figure  9a  -  Nonorthogonal 


Figure  9b  -  Orthogonal 


Figure  9  -  Coordinate  System  for  Region  with  Concave  Boundary 


the  orthogonal  mesh  must  have  rather  fine  spacing  near  the  concave  upper  boundary 
to  accommodate  the  curvature.  The  fine  mesh  spacing  raises  a  question  concerning 
the  possibility  of  a  singularity  in  this  vicinity.  This  is  a  valid  question,  one 
that  needs  to  be  considered  each  time  the  boundary-fitted  coordinate  technique  is 
used  since  almost  any  generating  system  can  produce  unacceptable  meshes  for  partic 
ular  regions.  To  verify  that  the  fine  mesh  spacing  in  Figure  9b  does  not  indicate 
a  singularity  in  the  transformation,  we  have  refined  the  mesh.  Figure  10  compares 
two  different  grids,  one  coarse  with  1681  points  and  the  other  fine  with  6561 
points,  generated  for  the  concave  region.  The  fact  that  corresponding  grid  lines 
are  in  about  the  same  position  in  both  meshes  confirms  that  the  coarse  discretiza¬ 
tion  yields  a  good  approximate  solution  to  the  exact  problem.  A  further  confirma¬ 
tion  comes  from  consideration  of  the  Jacobian  at  the  midpoint  of  the  upper  boun¬ 
dary.  The  value  of  J=0. 00186  computed  on  the  coarse  mesh  is  seen  to  agree  very 
well  with  J=0. 00047  on  the  fine  mesh  when  it  is  taken  into  account  that  these  quan 
tities  should  differ  by  a  factor  of  four  because  of  the  discretization  details. 
There  is  no  indication  of  a  zero  Jacobian  in  the  region. 

In  order  to  demonstrate  some  of  the  problems  that  can  arise,  we  attempted  to 
generate  an  orthogonal  mesh  on  a  region  similar  to  the  previous  one  but  with 
greater  curvature  of  the  concave  boundary.  The  grid  shown  in  Figure  11a  was  gen¬ 
erated  by  Equation  (2)  with  P  =  Q  =  0  while  the  unacceptable  grid  in  Figure  lib  was 
generated  by  the  system  of  Equation  (12).  Of  course  an  unacceptable  mesh  such  as 
this  one  with  crossing  lines  indicates  a  singular  transformation  which  can  often 
lead  to  numerical  difficulties.  But  problems  like  this  can  also  arise  from  itera¬ 
tive  schemes  based  on  the  Poisson  system  if  the  forcing  functions  are  not  chosen 
carefully.  To  verify  this  we  computed  directly  forcing  functions  P  and  Q  using 
Equation  (2)  and  x  and  y  given  as  in  Figure  lib.  We  then  solved  Equation  (2) 


Figure  10  -  Comparison  of  Orthogonal  Grids  for  Region  with  Concave  Bounds 

1681  Points,  —  —  -  6561  Points 


Figure  11  -  Concave  Region 


iteratively  for  x  and  y  using  this  P  and  Q,  thus  regenerating  the  grid  of  Figure 
lib. 


As  the  next  example,  consider  a  doubly-connected  region  bounded  by  concentric 
circles  as  shown  in  Figures  12  and  13.  Since  this  region  is  symmetric  to  a  line 
through  the  center,  each  grid  was  generated  for  half  the  region  and  reflected  in 
the  line  of  symmetry.  The  symmetry  line  was  treated  as  a  boundary  with  fixed  coor¬ 
dinate  distribution,  thus  assuring  that  $  =  0  at  the  corners  of  the  computational 
region.  The  spacing  on  the  outer  boundary,  but  not  on  the  inner  boundary,  was  uni¬ 
form.  Had  the  spacing  on  both  boundaries  been  uniform,  the  grid  produced  by  the 
Laplace  generating  system  (Figure  12a)  would  have  been  the  usual  polar  coordinate 
system  which  is  orthogonal.  In  Figure  12b,  the  line  of  symmetry  was  taken  as  a 
horizontal  line  through  the  center  of  the  figure  while  the  line  of  symmetry  for 
Figure  13  was  a  vertical  line  through  the  center.  Interestingly,  the  two  orthogo¬ 
nal  grids  thus  generated  (Figures  12b  and  13)  are  quite  dissimilar  as  a  result  of 
different  points  being  held  constant  after  the  same  initial  guess. 

For  a  final  test  of  the  method,  we  consider  a  square  region  with  the  coordi¬ 
nates  distributed  uniformly  on  three  boundaries,  nonuniformly  on  the  fourth.  Fig¬ 
ure  14a  shows  a  nonorthogonal  grid  with  441  points  and  Figure  14b  a  converged  solu¬ 
tion  to  Equation  (12)  with  3  ^  0.00020  throughout  the  field  except  on  the  boun¬ 
daries.  In  an  attempt  to  improve  the  orthogonality,  the  number  of  grid  points  was 
increased  to  1681  as  seen  in  Figure  15.  This  refinement  of  the  grid  did  not 
improve  the  orthogonality.  Thus  we  believe  that  no  orthogonal  mesh  exists  for  this 
configuration  and  that  the  existence  of  a  numerical  solution  with  a  small  non¬ 
zero  6  results  from  the  weak  connection  between  the  interior  grid  points  and  the 
corner  points,  the  only  points  where  8  *  0  is  enforced. 
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Figure  13  -  Orthogonal  Coordinate  System  for  Annular  Region 
(Vertical  Symmetry  Line) 
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Figure  15  -  Coordinate  System  with  Unequal  Spacing  on  One  Boundary  (1681  Points) 


SUMMARY 


Orthogonal  boundary-fitted  coordinate  systems  can  be  generated  for  a  specified 
boundary  coordinate  distribution  by  the  system  of  Equation  (2)  if  P  and  Q  can  be 
found  such  that  3=0.  One  way  to  do  this  would  be  to  solve  the  combined  system 
of  Equations  (2)  and  (6).  Iterative  solution  of  this  system  using  Equations  (9) 
diverged  on  a  simple  test  problem.  An  alternative  iterative  scheme  given  by  Equa¬ 
tions  (2)  and  (11)  converged  extremely  slowly  for  simple  test  problems.  An  itera¬ 
tive  method  based  on  the  solution  of  Equation  (12)  rather  than  on  the  Poisson  gen¬ 
erating  system  had  better  convergence  properties  and  was  used  to  generate  orthogo¬ 
nal  boundary-fitted  systems  for  several  test  problems. 

While  questions  remain  concerning  the  existence  and  uniqueness  of  orthogonal 
coordinate  systems,  and  although  there  is  much  room  for  improvement  in  the  methods 
for  developing  them,  the  generating  methods  presented  here  add  to  the  available, 
useful  techniques  for  constructing  these  systems.  Of  course,  it  must  be  realized 
that  a  particular  orthogonal  coordinate  system  is  not  necessarily  better  than  any 
given  nonorthogonal  system.  In  general,  it  seems  beneficial  to  require  orthogonal¬ 
ity  whenever  it  does  not  lead  to  a  serious  loss  of  other  desirable  coordinate  sys¬ 
tem  properties. 
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